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-The  development  of  a  panel  method  suitable  for  the  analysis  of  ducted  propellers  is  presented. 
The  method  is  first  applied  to  the  problem  of  the  two-dimensional  hydrofoil,  the  propeller  with  hub 
and  the  axisymmetric  duct  in  uniform  flow.  Some  comparisons  are  made  with  exact  solutions  and 
with  other  panel  codes.  The  difficulties  associated  with  the  modeling  of  ducted  propeller  flows  are 
discussed.  Convergence  of  the  method  for  a  typical  ducted  propeller  is  shown.  The  results  include 
overall  forces  on  the  propeller  and  the  duct,  circulation  distributions,  and  chordwise  pressure 
distributions  on  the  duct  at  different  positions  between  blades/  > 
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1.  Introduction 

The  design  of  ducted  propellers  contains  all  the  diffi¬ 
culties  associated  with  the  design  of  open  propellers,  in¬ 
cluding  the  conflicting  considerations  of  efficiency, 
cavitation,  strength,  and  vibration.  However,  one  now  has 
two  complex  geometries  to  deal  with,  thus  increasing  the 
number  of  variables  in  the  problem,  and  increasing  the 
complexity  of  the  hydrodynamic  problems  which  must  be 
dealt  with. 

As  with  the  case  of  open  propellers,  the  approaches 
available  to  the  designer  vary  greatly  in  complexity.  One 
can  simply  select  an  existing  geometry  from  published 
systematic  series  data,  or  one  can  develop  a  design  which 
meets  the  unique  requirements  of  a  particular  application. 
In  the  latter  case,  one  must  either  trust  in  the  accuracy 
of  an  existing  analytical  design  method,  conduct  extensive 
model  tests,  or  settle  on  some  compromise  between  the 
two. 

The  disadvantage  of  model  testing  is  that  it  is  both  time- 
consuming  and  expensive,  thus  limiting  the  extent  to 
which  a  wide  variety  of  design  alternatives  can  be  ex¬ 
plored.  This  is  certainly  true  for  open  propellers,  but  even 
more  true  for  ducted  propellers,  where  both  the  number 
of  options  to  be  studied  and  the  complexity  of  the  model 
test  is  much  greater.  Moreover,  the  ever-present  question 
of  scale  effect  errors  is  greater  for  ducted  propeller  tests 
due  to  the  lower  Reynolds  number  of  the  duct  as  compared 
to  that  of  a  rotating  propeller  blade. 

A  reliable  analytical  method  for  the  hydrodynamic  de¬ 
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sign  and  analysis  of  ducted  propellers  would  therefore  be 
desirable.  This  need  has  been  recognized  for  many  years, 
and  substantial  progress  toward  this  goal  has  been  made 
as  a  result  of  the  contributions  of  a  large  number  of  re¬ 
searchers  working  in  several  distinct  fields. 

The  largest  element  of  previous  research  in  ducted  pro¬ 
pellers  has  focused  on  the  axisymmetric  problem  of  a  duct 
interacting  with  actuator  disk,  or  equivalently,  with  an 
infinite-bladed  lifting  line.  The  earliest  representations  of 
the  duct  were  fully  linearized,  with  singularities  repre¬ 
senting  the  duct  thickness  and  loading  projected  onto  a 
mean  cylindrical  surface.  This  was  extended  to  place  the 
singularities  on  geometrically  more  complicated  surfaces, 
to  introduce  some  nonlinear  aspects  in  the  computation 
of  pressures,  and  finally,  to  represent  the  duct  more  exactly 
by  singularities  placed  on  the  actual  duct  surface.  The  work 
of  Gibson  and  Lewis  [l]4  and  of  Glover  and  Ryan  [2]  is 
representative  of  this  stage  of  ducted  propeller  theory. 
More  recently,  Falcao  de  Campos  [3]  again  retained  an 
actuator  disk  model  of  the  propeller,  but  introduced  a 
numerical  method  to  account  for  shear  flow  effects  caused 
by  a  radially  varying  inflow  field. 

At  the  same  time,  propeller  theory  was  advancing  from 
lifting  line  theory  to  well-developed  lifting-surface  theo¬ 
ries  which  could  account  for  the  increasingly  complex  ge¬ 
ometries  which  were  being  adopted  by  propeller 
designers.  Nevertheless,  relatively  few  attempts  to  com¬ 
bine  more  exact  representations  of  blades  and  ducts  have 
appeared  in  the  literature.  Two  recent  publications  which 
address  this  problem  are  by  VanHouten  [4]  and  Feng  and 
Dong  [5], 

Another  important  source  of  research  applicable  to 
ducted  propellers  comes  from  the  field  of  turbomachinery. 
In  this  case,  of  course,  the  flow  is  limited  to  the  interior 
of  a  prescribed  axisymmetric  body  where  mass  flow  con¬ 
tinuity  immediately  relates  the  upstream  and  downstream 
flows.  Hence  the  similarity  between  the  two  problems  is 


4  Numbers  in  brackets  designate  References  at  end  of  paper. 
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less  when  the  duct  is  short  relative  to  the  propeller  radius. 
On  the  other  hand,  the  local  flow  between  the  propeller 
tip  and  the  duct  must  be  fundamentally  very  similar.  Sam¬ 
ples  of  the  extensive  literature  in  this  field  are  papers  by 
Lakshminarayana  [6]  and  by  Booth,  Dodge  and  Hepworth 

m. 

The  development  of  methods  for  treating  axisymmetric 
shear  flows  is  also  an  important  aspect  of  turbomachinery 
theory,  and  much  of  this  is  applicable  to  the  ducted  pro¬ 
peller  problem.  An  example  is  the  streamline  curvature 
method  presented  by  McBride  [8], 

Finally,  the  broad  field  of  discretized  boundary  integral 
methods  has  application  to  the  specific  problem  of  the 
ducted  propeller. 

In  this  paper  we  concentrate  on  the  formulation  and 
numerical  solution  of  the  steady,  potential  flow  problem 
of  a  duct,  hub  and  propeller  with  a  finite  number  of  blades. 
Except  for  the  tip  flow  problem,  which  we  address  briefly, 
real  fluid  effects  will  not  be  considered.  Our  objective  is 
to  develop  a  robust  potential  flow  foundation  for  the  flow 
produced  by  geometrically  complex  ducted  propellers. 


2.  General  discussion  of  panel  methods 

Discretized  boundary  integral  methods,  more  com¬ 
monly  referred  to  as  panel  methods,  are  rapidly  becoming 
an  essential  tool  for  the  analysis  of  flows  around  complex 
three  dimensional  objects.  Since  the  publication  of  the 
pioneering  work  of  Hess  and  Smith  [9]  in  1964,  a  large 
number  of  different  panel  methods  have  appeared.  Some 
of  these  can  be  categorized  as  production  codes,  which 
can  be  used  for  a  wide  variety  of  applications.  Surveys  of 
some  of  the  major  production  codes  may  be  found  in  pub¬ 
lications  by  Hess  [10],  Hunt  [11],  Margason,  et  al  [12], 
Maskew  [13],  and  Youngren  et  al  [14].  Others  are  special 
purpose  codes  which  have  been  developed  with  specific 
applications  in  mind.  Examples  of  the  latter  are  the  pro¬ 
peller  duct  code  developed  by  Gibson  and  Lewis  [1]  which 
we  discussed  earlier,  and  a  recently  developed  propeller 
panel  code  presented  by  Hess  and  Valarezo  [15]. 

The  governing  equation  for  most  of  the  original  panel 
methods  was  the  requirement  that  the  normal  velocity  be 
zero  at  a  selected  point  on  each  panel.  Subsequently,  Mor- 
ino  [16]  introduced  a  panel  method  based  on  Green’s  for¬ 
mula  in  which  the  primary  unknown  was  the  potential. 
Both  velocity  and  potential  methods  are  in  current  use, 
and  each  has  its  advantages  and  disadvantages. 

In  addition  to  this  fundamental  difference,  various  com¬ 
binations  of  source,  vortex  and  dipole  singularities  are  em¬ 
ployed  in  the  various  methods.  This  might  seem  arbitrary, 
since  there  is  only  one  unique  solution  for  a  given  body 
geometry  and  inflow  field.  The  explanation  is  that  differ¬ 
ent  formulations  yield  different  flows  internal  to  the  body, 
but  identical  flows  in  the  exterior  field.  Since  the  internal 
flow  has  no  physical  significance,  different  distributions  of 
singularities  on  the  boundaries  of  a  body  can  produce  same 
correct  physical  flow.  A  brief  but  rigorous  derivation  of 
the  equivalence  of  velocity  and  potential  formulations,  and 
of  the  relationship  between  the  types  of  boundary  singu¬ 
larities  is  presented  in  Appendix  1.  This  can  be  summa¬ 
rized  with  the  happy  conclusion  that  everybody  is  right! 

A  numerical  implementation  of  any  of  these  fondamen- 
tal  approaches  involves  approximations  of  various  sorts, 
and  it  is  here  that  various  methods  differ.  The  body  surface 
can  be  approximated  by  plane  quadrilateral  elements  with 
constant  singularity  density  within  each  panel,  or  by 
curved  panels  with  varying  singularity  densities.  Methods 


employing  the  former  are  termed  low  order  panel  methods 
while  the  latter  are  referred  to  as  high  order  methods. 
High  order  methods  can  achieve  a  prescribed  level  of 
accuracy  with  fewer  panels,  but  the  code  is  more  com¬ 
plicated  and  the  required  amount  of  computation  per 
panel  is  higher.  Since  the  relative  merits  of  low  and  high 
order  methods  depend  both  on  the  specific  problem  and 
on  the  fundamental  method  being  used,  it  is  not  surprising 
that  the  issue  is  controversial. 

Finally,  problems  involving  lift  require  the  imposition 
of  a  Kutta  condition  at  the  trailing  edge.  Implementation 
of  the  Kutta  condition  in  a  discretized  problem  again  in¬ 
volves  approximations,  and  these  vary  between  methods. 

A  comprehensive  study  of  the  field  of  panel  methods 
has  recently  been  conducted  by  J.-T.  Lee  [17]  in  order  to 
determine  the  most  suitable  formulation  for  application 
to  marine  propeller  and  duct  problems.  The  conclusion 
was  that  a  low  order  potential  based  panel  method,  cou¬ 
pled  with  nonlinear  pressure  Kutta  condition,  would  be 
best.  Some  of  the  principal  reasons  for  this  conclusion  are, 

•  While  practically  all  methods  work  well  for  thick  sec¬ 
tions,  the  potential  method  is  substantially  more  ac¬ 
curate  for  very  thin  sections.  This  is  particularly 
important  for  marine  propellers,  where  thickness/ 
chord  ratios  typically  vary  from  20  percent  at  the  root 
to  as  little  as  2  percent  near  the  tip. 

•  The  influence  coefficients  for  a  potential  induced  by 
unit  source  and  dipole  distributions  are  one  order  less 
singular  than  the  corresponding  influence  coefficients 
for  velocity.  As  a  result,  potential  based  methods  are 
expected  to  be  less  sensitive  to  errors  caused  by  ir¬ 
regular  paneling. 

•  The  computation  of  the  panel  influence  functions, 
which  is  a  major  contributor  to  the  total  computing 
effort,  is  faster  for  a  potential  method  than  for  a  ve¬ 
locity  method.  Of  course,  one  must  eventually  deter¬ 
mine  velocities  in  order  to  obtain  the  pressure 
distribution.  However,  this  can  be  accomplished  in  a 
more  accurate  and  efficient  way  by  differentiating  the 
potential  than  by  direct  computation  of  the  velocities 
induced  by  the  panel  singularities. 

•  Since  the  potential  influence  coefficients  are  scalar 
quantities,  the  total  storage  required  is  one  third  as 
great  as  for  a  velocity  method. 

•  A  low  order  potential  method  is  substantially  more 
accurate  than  a  low  order  velocity  method  for  the 
computation  of  internal  flows.  This  problem  is,  of 
course,  of  direct  concern  in  the  ducted  propeller  prob¬ 
lem. 

•  A  low  order  potential  method  is  probably  computa¬ 
tionally  more  efficient  than  a  high  order  velocity  or 
potential  method  for  incompressible  flows.  Even  if  this 
were  to  be  proved  subsequently  not  to  be  the  case, 
we  would  never  know  unless  we  had  an  efficient  low 
order  code  to  serve  as  a  basis  for  comparison. 

A  panel  calculation  for  the  flow  around  a  body  with  lift 
involves  the  fol’ owing  relatively  independent  steps, 

•  grid  generation, 

•  development  of  logic  to  take  advantage  of  symmetry, 

•  computation  of  the  panel  singularity  influence  func¬ 
tions, 

•  solution  of  the  simultaneous  equations  for  the  un¬ 
known  singularity  strengths, 

•  determination  of  local  velocities  and  pressures,  and 

•  determination  of  total  forces  and  moments. 

While  the  major  emphasis  in  this  paper  is  on  the  ducted 
propeller  problem,  this  represents  systematic  evolution  of 
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the  application  of  panel  methods  to  two-dimensional  hy¬ 
drofoils,  propellers  with  hubs,  isolated  ducts  and,  finally, 
propellers  with  hubs  and  ducts.  We  will  therefore  describe 
how  each  of  these  steps  is  accomplished  as  the  present 
method  is  applied  to  all  four  of  these  problems.  While  we 
will  not  be  discussing  the  structure  of  the  computer  codes, 
it  is  important  to  note  that  all  of  these  problems  contain 
common  elements  which  can  be  individually  programmed 
and  combined  in  a  way  which  is  most  efficient  for  the 
problem  at  hand. 


3.  The  two-dimensional  hydrofoil 

The  essential  first  step  in  a  panel  procedure  is  the  gen¬ 
eration  of  accurate  coordinates  for  the  vertices  of  each  of 
the  panels  representing  the  body.  Any  unfairness  intro¬ 
duced  by  an  imperfect  interpolation  process  will  result  in 
a  bumpy  pressure  distribution.  Most  hydrofoil  sections  are 
defined  by  the  superposition  of  a  mean  line  and  a  thickness 
form,  and  as  long  as  both  are  defined  analytically,  or  are 
approximated  by  cubic  splines  passing  through  an  accurate 
set  of  base  points,  the  required  degree  of  smoothness  is 
easy  to  achieve.  Sections  generated  by  conformal  mapping 
are  frequently  useful  to  compare  numerical  and  exact  so¬ 
lutions.  In  this  case,  the  evaluation  of  the  section  at  the 
desired  grid  points  can  be  carried  out  to  any  desired  de¬ 
gree  of  precision. 

The  next  step  is  to  devise  a  systematic  way  of  deter¬ 
mining  the  spacing  of  the  panels,  given  their  total  number. 
For  a  two-dimensional  hydrofoil,  an  arrangement  com¬ 
monly  referred  to  as  cosine  spacing  is  a  logical  choice.  For 
a  given  total  number  of  panels,  N,  the  mean  line  ordinate 
and  thickness  is  first  evaluated  at  the  following  points  along 
the  nose-tail  line, 

'-§(■—  (^))' 

j  =  1,  2, ... ,  NI2  +  I  (1) 

where  x  is  the  coordinate  along  the  nose-tail  line  of  a 
section  of  chord  length  c.  The  panel  boundaries  are  then 
obtained  by  adding  and  subtracting  the  half-thickness  of 
the  section  at  right  angles  to  the  mean  line.  This  concen¬ 
trates  the  elements  at  the  leading  and  trailing  edges, 
where  greater  resolution  is  required.  This  characteristic 
can  be  seen  from  the  example  shown  in  Fig.  1. 

The  general  two-dimensional  hydrofoil  problem  does 
not  possess  any  symmetry  which  can  be  used  to  increase 
the  efficiency  of  the  computation.  This  would  not  be  an 
important  issue  in  any  case,  since  two-dimensional  calcu¬ 
lations  use  such  a  small  amount  of  computer  time. 

The  exact  influence  functions  for  two-dimensional  con¬ 
stant  strength  source  and  dipole  elements  are  well  known 
[18].  Unlike  three-dimensional  applications,  where  sub- 


Flfl.1  awdwtoedtotrtoutlon  of  panels  produced  by  oosfoeapecing 
as  defined  in  equation  (1).  The  two-dimensional  section  hat  a  thick¬ 
ness  chord  ratio  of  5  percent  The  number  of  panels  is  40 


stantial  savings  in  computing  time  can  be  achieved  by  the 
use  of  far  field  approximations  for  distant  panels,  the  time 
savings  for  a  two-dimensional  calculation  are  not  worth 
the  trouble. 

An  approximate  Kutta  condition  for  a  potential  flow 
panel  solution,  first  introduced  by  Morino  [16],  required 
that  the  strength  of  the  dipole  sheet  in  the  wake  be  equal 
to  the  difference  in  the  value  of  the  dipole  strength  of  the 
two  panels  adjacent  to  the  trailing  edge.  This,  combined 
with  the  discretized  statement  of  Green’s  formula  given 
in  Appendix  1,  results  in  a  system  of  linear  equations  for 
the  unknown  dipole  strengths  of  each  panel.  However,  J.- 
T.  Lee  [17]  found  that  this  form  of  the  Kutta  condition 
contained  a  fundamental  error  when  the  free  stream  con¬ 
tained  a  component  in  the  direction  of  a  line  connecting 
the  control  points  of  the  two  trailing  edge  panels.  This 
finding  led  to  a  study  of  the  implementation  of  the  Kutta 
condition  for  a  potential  based  panel  method.  The  con¬ 
clusion  was  that  the  most  accurate  and  reliable  Kutta  con¬ 
dition  was  a  requirement  of  equal  pressures  on  the  two 
trailing  edge  elements.  This  is  similar  to  the  conclusion 
reached  by  Hess  [10]  for  his  velocity  method.  It  should  be 
noted,  however,  that  the  differences  between  the  results 
obtained  by  the  original  Morino  Kutta  condition  and  the 
present  pressure  Kutta  condition  are  very  small  for  thin 
two-dimensional  sections.  The  differences  become  signif¬ 
icant  for  thick  sections  and  for  three-dimensional  flows. 

This  introduces  a  nonlinear  aspect  to  the  solution.  How¬ 
ever,  an  iterative  scheme  can  be  employed  whereby  the 
initial  solution  is  obtained  using  the  original  Morino  con¬ 
dition,  and  successive  wake  dipole  strengths  are  adjusted 
based  on  the  error  in  the  computed  trailing  edge  element 
pressures. 

The  coefficient  matrix  is  unchanged  in  the  process.  For 
the  small  number  of  unknowns  involved  for  two-dimen¬ 
sional  flows,  a  Gauss  reduction  can  be  used  at  the  outset 
to  decompose  the  matrix  into  lower  and  upper  triangular 
forms  for  subsequent  solutions  with  different  right  hand 
sides.  The  converged  solution  for  the  dipole  strengths  then 
represents  the  distribution  of  potential  around  the  surface 
of  the  hydrofoil. 

Surface  velocities  may  then  be  obtained  either  by  nu¬ 
merical  differentiation  of  the  potential,  or  by  direct  cal¬ 
culation  of  the  source  and  dipole  panel  velocity  influence 
functions.  The  latter  approach  is  generally  not  as  accurate, 
since  the  velocity  influence  functions  are  more  singular 
and  therefore  more  sensitive  to  the  position  of  the  control 
points  within  each  panel.  This  problem  is  similar  to  that 
of  a  vortex  lattice  calculation,  where  the  position  of  the 
control  points  is  known  to  be  critical.  As  a  result,  numerical 
differentiation  of  the  potential,  with  all  its  inherent  dan¬ 
gers,  is  still  the  better  alternative. 

The  crudest  form  of  numerical  differentiation  simply 
consists  of  representing  the  velocity  at  the  panel  bound¬ 
aries  as  the  difference  in  potential  of  the  two  adjacent 
elements  divided  by  the  distance  between  their  centers. 
A  more  accurate  scheme  is  to  employ  either  a  cubic  spline 
interpolation  procedure,  or  a  second-order  finite  differ¬ 
ence  method,  with  arc  length  along  the  surface  as  a  pa¬ 
rameter.  We  have  chosen  this  last  approach,  although  we 
have  succeeded  in  obtaining  convergent  results  with  both 
of  the  alternatives  just  discussed. 

The  final  step  is  the  computation  of  total  forces.  We 
know,  of  course,  that  for  two-dimensional  in  viscid  flow  the 
drag  is  zero  and  the  lift,  from  Kutta-Joukowski's  law,  is 
pUr.  Moreover,  the  circulation,  T,  is  equal  to  the  known 
potential  jump  or  dipole  strength  in  the  wake.  However, 
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Flfl.  2  Chordwise  distribution  of  pressure  coefficient  C„  for  a  2  percent  thick 
symmetrical  Karman-Trefftz  section  at  an  angle  of  attack  of  2  deg.  The  trailing 
edge  angle  specified  in  the  mapping  function  is  3.6  deg.  Results  obtained  with 
the  present  panel  method  are  given  for  20.  40  and  SO  panels 


to  practice  for  more  complicated  three-dimensional  flows, 
we  will  also  calculate  the  force  the  hard  way  by  integrating 
the  pressure  around  the  section.  This  can  be  implemented 
by  multiplying  the  pressure  computed  at  the  midpoint  of 
each  panel  by  the  panel  arc  length,  forming  a  set  of  N 
concentrated  force  vectors  in  the  direction  of  the  surface 
normals,  and  adding  them  up.  The  error  incurred  by  local 
pressure  integration  can  then  be  readily  determined. 

A  typical  pressure  distribution  for  a  thin  section  is  given 
in  Fig.  2.  Results  are  given  for  a  Karman-Trefftz  section 
whose  exact  solution  can  be  obtained  by  conformal  map¬ 
ping  [19].  The  numerical  results  are  all  very  close  to  the 
exact  solution,  even  with  20  panels.  However,  40  panels 
would  probably  provide  a  better  definition  of  the  pressure 
near  the  leading  edge.  Results  for  thicker  sections  will  not 
be  shown  here  since  the  agreement  between  exact  and 
numerical  results  is  even  better. 

Table  1  provides  a  comparison  of  the  computed  values 
of  lift  and  drag  coefficient  obtained  by  summation  of  panel 
forces  with  exact  values  obtained  by  conformal  mapping. 
This  is  a  demanding  test  since  the  thin  sections  develop 


Table  1  Effect  of  the  number  of  chordwise  pends  on  com¬ 
puted  HR  and  drag  coefficient  The  exact  values  are  Ct  =  0.223 
C„  =  0 


N 

Co 

Cl 

20 

0  0061 

0  234 

40 

00043 

0  329 

80 

00030 

0  226 

very  sharp  pressure  peaks  at  the  leading  edge.  Using  80 
elements,  the  error  in  drag  is  less  than  1  percent  of  the 
lift.  A  high  level  of  accuracy  can  therefore  be  achieved, 
but  one  must  be  willing  to  use  a  large  number  of  panels. 

4.  The  propeller  and  hub 

The  problem  of  grid  generation  for  a  propeller  blade 
can  be  separated  into  two  parts — chordwise  spacing  anc’ 
radial  spacing.  The  former  is  identical  to  the  two-dimen¬ 
sional  problem,  and  we  have  again  adopted  cosine  spacing, 
as  illustrated  in  Figs.  3  and  4. 

The  best  choice  of  the  radial  distribution  of  panel  size 
for  a  propeller  blade  is  not  obvious.  For  the  example  shown 
we  have  used  a  spacing  algorithm  which  concentrates  the 
elements  at  the  tip.  In  this  case,  the  radii  for  M  panels  are, 

rm  =  Rh  +  (R  -  Rh)sin  ~  1>), 

m  =  I,  2 . M  +  1  (2) 

where  Rh  is  the  hub  radius,  R  is  the  tip  radius,  and  M  is 
the  number  of  panels  over  the  radius.  For  some  applica¬ 
tions,  however,  greater  resolution  near  the  hub  might  be 
desirable,  in  which  case  cosine  spacing  rather  than  the 
spacing  given  in  equation  (2)  might  be  preferred. 

The  propeller  hub  and  its  intersection  with  the  propeller 
blades  is  more  complicated.  However,  we  have  assumed 
for  the  present  that  the  radius  of  the  hub  is  constant  be¬ 
tween  the  leading  and  trailing  edges  of  the  blade,  and  can 
be  approximated  by  a  quartic  function  of  axial  distance 
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Fig.  3  Panel  arrangement  viewed  from  upstream  for  a  three-bladed 
propeller  with  a  hub.  The  propeller  is  DTNSRDC-41 18,  whose  char¬ 
acteristics  may  be  found  in  [23] 

downstream.  The  upstream  portion  of  the  hub  is  treated 
as  a  semi-infinite  cylinder.  Hub  definition,  and  the  hub- 
blade  intersection  problem,  is  then  made  much  easier. 

The  axial  spacing  of  the  panels  on  the  hub  matches  that 
of  the  blade  in  the  interval  between  the  leading  and  trail¬ 
ing  edges,  and  the  circumferential  interval  along  the  hub 
is  divided  into  a  prescribed  set  of  equal  intervals.  This 
generates  a  more  or  less  helical  pattern  of  paneling  on  the 
hub.  The  axial  coordinates  of  the  intermediate  panels  must 
be  adjusted,  particularly  near  the  leading  edge,  in  order 
to  avoid  badly  shaped  panels  in  this  region.  If  the  axial 
coordinates  were  all  required  to  be  the  same,  panels  near 
the  leading  edge  could  actually  turn  inside-out.  This  ar¬ 
rangement  produced  some  interesting  hub  pressure  dis¬ 
tributions  which  caused  us  to  take  a  closer  look  at  our  hub 
panel  graphics.  The  result  was  the  addition  of  an  automatic 
paneling  scheme  which  appears  to  produce  reasonable 
panel  shapes. 

The  paneling  on  the  hub  upstream  of  the  propeller  lead¬ 
ing  edge  is  purely  helical,  with  a  pitch  matching  the  root 
section  pitch  of  the  propeller.  The  axial  spacing  is  half 
cosine,  with  fine  spacing  near  the  propeller  and  coarse 
spacing  upstream.  The  arrangement  on  the  hub  down¬ 
stream  of  the  trailing  edge  is  similar,  except  that  the  pitch 
is  required  to  match  the  corresponding  pitch  of  the  trailing 
vortex  wake  at  the  hub  intersection. 

As  with  any  three-dimensional  lifting  surface,  a  trailing 
vortex  wake  extends  downstream  from  the  blades  and  hub, 
which  must  also  be  paneled  with  quadrilateral  dipole  ele¬ 
ments.  We  will  use  a  model  for  the  trailing  vortex  wake 
identical  to  that  developed  by  Greeley  and  Kerwin  [20] 
for  use  with  a  vortex  lattice  propeller  analysis  code.  Stated 
briefly,  the  model  divides  the  wake  into  a  transition  wake 
extending  from  the  trailing  edge  to  a  specified  distance 
downstream.  In  this  region,  the  wake  is  represented  as  a 
discretized  vortex  sheet  with  a  prescribed  contraction 
shape  and  an  axially  varying  pitch.  The  second  region, 
designated  the  ultimate  wake,  is  modeled  more  crudely 
by  a  single  concentrated  tip  vortex  from  each  blade,  and 
a  hub  vortex.  The  quadrilateral  dipole  panels  in  the  wake 


are  geometrically  identical  to  the  vortex  lattice  elements 
used  in  [201 

We  turn  next  to  the  question  of  symmetry.  For  the 
steady  flow  problem,  both  the  geometry  and  the  loading 
is  repeated  identically  on  each  blade  and  on  each  inter¬ 
blade  segment  of  the  hub.  The  number  of  unknown  dipole 
strengths  to  be  determined  is  therefore  the  total  number 
of  elements  (as  seen,  for  example,  in  Fig.  3)  divided  by  the 
number  of  blades.  In  fact,  the  geometry  is  only  determined 
for  one  blade,  and  the  influence  functions  incorporate  a 
built-in  summation  of  the  influence  of  a  set  of  source  and 
dipole  panels  with  angular  coordinates  incremented  by 
the  angle  between  blades. 

This  is  equivalent  to  having  equal  paneling  on  all  blades. 
Our  earlier  vortex  lattice  codes  [20,  21]  used  coarser  spac¬ 
ing  on  the  other  blades,  together  with  an  interpolation 
scheme  to  account  for  their  influence  on  the  first,  or  key, 
blade.  However,  the  panel  influence  coefficient  algorithm 
which  we  are  using  contains  an  efficient  far  field  approx¬ 
imation  which  reduces  the  time  required  to  compute  dis¬ 
tant  elements.  Therefore  the  substantial  coding 
complications  required  to  handle  different  spacing  on 
other  blades  was  not  considered  to  be  necessary  in  this 
case. 

The  influence  functions  for  quadrilateral  source  and  di¬ 
pole  panels  are  computed  using  the  formulation  developed 
by  Newman  [22].  The  potential  is  computed  exactly  for 
nearby  panels,  while  it  is  approximated  by  a  multipole 
expansion  for  more  distant  panels.  Finally,  the  panels 
which  are  sufficiently  distant  are  treated  as  point  sources 
and  dipoles.  This  is  really  the  heart  of  any  panel  code,  and 
it  is  essential  that  the  computer  code  for  the  influence 
function  be  both  robust  and  efficient. 

The  Kutta  condition  employed  in  the  propeller  problem 
is  logically  similar  to  that  used  for  the  two-dimensional 
hydrofoil.  We  found  that  an  explicit  pressure  Kutta  con¬ 
dition  was  particularly  necessary  near  the  tip  where  large 
radial  induced  velocities  combined  with  locally  swept  trail¬ 
ing  edges  combined  to  produce  inaccurate  results.  A  de¬ 
tailed  development  of  the  Kutta  condition  may  be  found 
in  Appendix  2. 

The  resulting  system  of  simultaneous  equations  is  similar 
to  that  described  for  a  simple  two-dimensional  hydrofoil, 
but  the  total  number  of  unknowns  is  much  greater.  As  a 
result,  the  time  required  for  a  Gauss  elimination  would 
start  to  become  prohibitive,  and  an  iterative  matrix  solver 


Fig.  4  Panel  arrangement  viewed  from  downstream  for  a  three- 
bladed  propeller  with  a  hub.  The  propeller  is  DTNSRDC-41 1 8,  whose 
characteristics  may  be  found  in  (23) 
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is  more  efficient.  We  have  employed  an  accelerated  it¬ 
erative  matrix  solver  following  the  procedure  developed 
by  Clark  [23].  While  this  was  developed  for  use  with  Hess’s 
velocity  based  panel  codes,  we  have  found  that  it  con¬ 
verges  very  rapidly  for  the  kind  of  coefficient  matrix  en¬ 
countered  with  a  potential  method.  For  example,  the 
solution  for  2000  unknowns,  using  a  DEC  Micro vax  II, 
takes  fifteen  minutes,  including  the  time  required  for  ex¬ 
ternal  storage  and  retrieval  of  the  matrix  elements.  On 
the  other  hand,  a  Gauss  reduction,  if  carried  out  entirely 
in  memory,  would  require  an  estimated  time  of  one  hour. 
Since  memory  limitations  would  make  the  latter  option 
impossible  for  such  a  large  matrix,  external  storage  would 
also  be  required,  and  the  resulting  time  would  be  much 
more  than  one  hour. 

The  calculation  of  pressure  distributions  is  again  carried 
out  by  determining  the  surface  velocity  by  differentiating 
the  potential.  However,  in  this  case  a  two-dimensional 
interpolation  is  required.  As  shown  in  Appendix  3,  this  is 
accomplished  by  determining  the  derivatives  of  the  po¬ 
tential  along  two  nonorthogonal  curvilinear  coordinates 
formed  by  the  blade  and  hub  paneling.  The  magnitude  of 
the  velocity  is  then  determined,  which  can  then  be  used 
to  determine  a  nondimensional  pressure  coefficient, 


Blade  pressure  distributions  at  three  radii  are  shown  in 
Fig.  5,  while  pressure  distributions  along  the  hub  are 
shown  in  Fig.  6. 


Finally,  the  total  propeller  thrust  and  torque  can  be 
obtained  by  summation  of  individual  panel  force  vectors. 
In  order  to  obtain  practically  useful  results  which  can  be 
compared  with  experiments,  viscous  drag  forces  are  added 
following  the  same  approximate  procedure  used  in  [20]. 
Tangential  force  vectors  are  added  in  each  panel  with 
magnitudes  equal  to  the  product  of  the  local  dynamic 
pressure,  the  panel  area  and  a  specified  viscous  drag  coef¬ 
ficient.  Thrust  and  torque  coefficients  obtained  in  this  way 
for  the  sample  propeller  at  several  advance  coefficients 
are  plotted  in  Fig.  7  together  with  experimental  results 
[24],  The  agreement  appears  to  be  satisfactory  in  this  case. 

5.  Axisymmetric  ducts  in  uniform  flow 

The  geometry  for  an  axisymmetric  duct  is  particularly 
simple,  since  one  only  needs  to  specify  a  meridional  section 
in  the  same  way  as  a  two-dimensional  hydrofoil.  The  chord- 
wise  distribution  of  panels  is  the  same  as  for  the  two- 
dimensional  case,  while  the  circumferential  paneling  is 
formed  from  constant  angular  increments. 

For  the  particular  case  of  uniform  inflow,  the  solution 
is  axisymmetric,  so  that  the  number  of  unknowns  is  equal 
to  the  number  of  chordwise  panels.  The  influence  func¬ 
tions  are  therefore  formed  by  summing  the  individual 
panel  contributions  circumferentially.  AD  other  steps  in 
the  solution  are  carried  out  in  the  same  way  as  for  a  two- 
dimensional  flow. 

The  first  example  shown  is  for  a  torus.  Fig.  8.  WhDe  this 
is  obviously  not  a  serious  candidate  for  a  propeUer  duct, 
an  analytical  solution  for  the  potential  was  developed  in 
1893  by  Dyson  [25].  The  analytical  solution  is  developed 


Fig.  S  Computed  chordwtse  pressure  distributions  for  DTNSRDC  propeller  4118 
operating  at  an  advance  coefficient  J  =  0.833.  The  pressure  coefficient  in  this 
case  is  normalized  on  the  uniform  advance  speed,  VA 
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Fig.  6  Computed  hub  pressure  distributions  for  DTNSRDC  propeller  4118  op¬ 
erating  at  an  advance  coefficient  J  =  0.833.  The  pressure  coefficient  in  this  case 
is  normalized  on  the  uniform  advance  speed,  VA.  The  pressure  distributions  are 
given  along  three  helical  panels  between  blades 


Fig.  7  Comparison  of  measured  and  calculated  thrust  and  torque 
coefficients  for  DTNSRDC  propeller  41 18.  The  experimental  results 
are  from  [23],  The  calculated  results  are  obtained  by  summation  of 
panel  pressures,  with  viscous  effects  treated  as  In  [19] 


Fig.  8  Panel  arrange  men*  for  toroidal  duct  with  20  circumferential 
and  20  chordwise  panels.  The  ratio  of  duct  cross  section  radius  to 
mean  radius,  rc/R  =  0.25 
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ng.  9  Distribution  of  potential  over  the  inner  and  outer  surfaces  of  a  toroidal  duct. 
The  exact  solution  due  to  Dyson  [24]  is  graphically  identical  to  the  80  x  80  panel 

solution 


as  an  expansion  in  powers  of  the  ratio  of  the  cross-sectional 
radius,  rc,  to  mean  radius,  R,  of  the  torus.  Since  the  circular 
cross  sections  of  the  duct  do  not  possess  a  natural  trailing 
edge,  the  circulation  is  made  unique  by  specifying  the 
angular  coordinate  of  the  downstream  stagnation  point. 

Figure  9  shows  the  computed  distribution  of  potential 
for  different  numbers  of  panels  for  a  torus  with  rJR  = 
0.25.  The  downstream  stagnation  point  is  specified  to  be 
at  an  angular  coordinate  of  zero  relative  to  the  axis  of  the 
torus. 

It  is  evident  that  the  computed  potential  distribution 
over  the  outer  surface  is  essentially  identical  for  the  three 
grids  tested.  The  result  for  the  inner  surface  is  not  con¬ 
verged  for  the  case  of  36  chordwise  and  9  circumferential 
panels,  but  the  results  for  the  two  finer  grids  are  almost 
identical.  The  analytical  solution  is  indistinguishable  from 
the  converged  numerical  results. 

The  next  case  is  a  very  long  duct  formed  from  an  NACA- 
0010  section  set  at  zero  angle  of  attack  on  a  mean  radius 
of  one  tenth  of  the  chord.  Results  for  this  duct  obtained 
with  several  panel  codes  were  given  by  Bristow  [26],  and 
additional  results  were  presented  by  Miranda  [27]  and  by 
Hess  [28]  in  a  discussion  to  Miranda’s  paper.  Predicting 
the  pressure  distribution  for  such  an  extreme  duct  is  a 
very  demanding  test  of  a  panel  code.  The  mass  flow 
through  the  duct,  which  is  imposed  by  the  Kutta  condition, 
is  extremely  high  in  this  case,  and  there  is  a  tendency  for 
all  panel  methods  to  underestimate  its  value.  This  is,  of 
course,  not  a  practical  issue  since  viscous  effects  would 
completely  invalidate  the  potential  flow  solution. 

Our  paneling  for  this  example  is  shown  in  Fig.  10,  and 
the  pressure  distributions  for  various  grids  are  shown  in 
Fig.  11.  The  first  four  calculations  are  for  a  fixed  number 
of  chordwise  panels  equal  to  36,  but  with  the  number  of 
circumferential  panels  varying  between  9  and  60.  The 
results  for  36  and  60  circumferential  panels  are  almost 
identical,  and  indicate  a  minimum  pressure  coefficient  of 
—  11.3.  Increasing  the  number  of  chordwise  elements  to 
60  reduces  the  minimum  pressure  coefficient  to  —12.5. 


Fig.  10  Panel  arrangement  for  nacelle  formed  from  an  NACA  0010 
section  with  a  chord/mean  radius  ratio  of  1 0.  There  are  36  chordwise 
and  18  circumferential  panels  in  this  illustration 


Some  results  obtained  with  other  panel  codes  are  shown 
in  Fig.  12,  which  is  taken  from  [27].  The  “exact”  solution, 
obtained  with  a  special  high  order  axisymmetric  panel 
code  [29]  shows  a  minimum  pressure  coefficient  of  —13.8 
which  is  slightly  lower  than  the  minimum  value  which  we 
obtained.  The  results  obtained  by  a  high  order  panel  code 
developed  by  Hess  [30],  and  by  the  QUADPAN  [27]  low 
order  potential  based  code,  are  similar  to  ours.  The  results 
obtained  by  the  Hess  low  order  velocity  based  code  [31] 
show  pressure  minima  much  closer  to  zero,  which  is  evi¬ 
dently  an  inherent  problem  with  that  method  for  internal 
flows. 

The  next  example  is  a  duct  with  proportions  which  are 
more  typical  ftr  a  ducted  propeller  application.  The  duct 
section  is  again  an  NACA-0015,  but  the  duct  radius  at  the 
leading  edge  is  now  equal  to  the  duct  chord.  The  duct 
sections  are  set  at  an  angle  of  attack  of  10  deg.  The  grid 
for  this  duct  is  shown  in  Fig.  13  and  pressure  distributions 
are  shown  in  Fig.  14  for  varying  numbers  of  panels.  Com¬ 
paring  this  with  Fig.  1 1  we  see  that  the  results  are  much 
more  convergent.  An  additional  measure  of  accuracy  is 
the  computed  drag  coefficient,  which  should  be  zero  in 
this  case.  Table  2  gives  the  values  of  drag  obtained  in  this 
case  for  different  numbers  of  panels.  The  computed  drag 
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Fly.  11  Pressure  distributions  for  nacelle  obtained  with  the  present  panel  code  with  different 

panel  arrangements 


becomes  smaller  as  the  number  of  panels  is  increased,  and 
it  is  somewhat  a  matter  of  judgment  to  set  an  acceptable 
value  for  the  drag  error.  If  this  duct  were  to  operate  with 
a  propeller,  the  drag  error  resulting  from  a  60  X  60  grid 
would  be  on  the  order  of  3  percent  of  the  duct  thrust,  or 
less  than  1  percent  of  the  total  duct  and  propeller  thrust. 


6.  Ducted  propellers 

We  now  turn,  finally,  to  the  problem  of  the  ducted 
propeller.  The  duct  and  hub  will  be  considered  logically 
as  one  unit  with  axisymmetric  geometry.  The  propeller 
blades,  however,  are  represented  by  a  vortex  lattice  using 
a  modification  of  the  procedure  described  in  [20].  This 
was  done  in  order  to  reduce  the  total  number  of  panels 
required  in  the  solution,  and  would  seem  justified  on  the 
basis  that  the  tip  sections  of  a  propeller  are  generally  much 
thinner  than  the  duct.  However,  we  plan  in  the  future  to 
provide  for  the  additional  alternative  of  representing  the 
blades  by  surface  panels  as  described  earlier. 


Fly.  12  Pressure  distributions  for  nscsNe  obtained  by  dHferent 
panel  methods  as  presented  by  Miranda  [26] 


The  procedure  works  as  follows:  The  flow  around  the 
combination  of  the  axisymmetric  hub  and  duct  is  first  com- 
puted,  using  exactly  the  same  method  as  for  the  duct  alone. 
The  velocity  induced  by  the  hub  and  duct  is  then  com¬ 
puted  as  a  set  of  field  points  in  the  interior  of  the  duct 
whose  coordinates  correspond  to  the  positions  of  the  con¬ 
trol  points  in  the  vortex  lattice  representation  of  the  pro¬ 
peller  blades.  The  strengths  of  the  vortex  lattice  elements 
on  the  propeller  blades  are  now  computed  for  this  partic¬ 
ular  spatially  varying  inflow  field  based  on  the  require¬ 
ment  that  the  total  normal  velocity  vanish  at  each  of  the 
control  points. 

The  perturbation  potential  induced  by  the  propeller 


Fly.  13  Panel  arrangement  for  a  duct  formed  from  an  NACA  0015 
section  with  a  chord /leading  edge  radkia  ratio  of  1.  The  duct  sections 
are  set  at  an  angle  of  attack  of  10  dag 
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Fig.  14  Pressure  distributions  for  duct  illustrated  in  Rg.  13  obtained  by  the 
present  panel  method  with  different  numbers  of  elements 


blades  can  now  be  computed  at  the  centroid  of  each  of 
the  panels  on  the  duct  and  hub.  This  modifies  the  right- 
hand  side  of  the  equations  determining  the  dipole 
strengths  on  the  hub  and  duct,  and  a  new  solution  for  the 
latter  is  obtained.  At  this  point  the  flow  around  the  duct 
and  hub  is  not  longer  axisymmetric. 

New  field  point  velocities  induced  by  the  duct  and  hub 
on  the  propeller  are  then  computed,  and  the  propeller 
solution  is  repeated.  The  iterations  continue  until  the 
changes  to  both  the  duct  and  hub  and  to  the  propeller  are 
within  a  prescribed  tolerance.  This  is  generally  accom¬ 
plished  within  five  or  six  iterations,  as  will  be  shown  later 
for  a  specific  example. 

Substantial  complications  arise  in  this  problem  as  a  result 
of  the  need  for  compatibility  between  the  grids  on  the 
duct  and  hub  and  on  the  blades.  For  example,  if  the  panels 
on  the  duct  lie  in  meridional  planes,  as  was  the  case  for 
the  axisymmetric  problem,  the  duct  grid  would  intersect 
the  tip  panels  of  the  propeller  blades  and  their  trailing 
vortex  wakes.  The  panels  on  the  duct  and  hub  are  there¬ 
fore  helical,  and  exactly  match  the  blade  and  trailing  vor¬ 
tex  grid  at  the  tip  and  at  the  hub. 

This  seemed,  superficially,  to  be  just  a  matter  of  ge¬ 
ometry  which  should  not  be  too  difficult  to  handle.  We 
there  tore  embarked  on  a  necessary  test  of  the  numerics, 
in  which  we  calculated  an  axisymmetric  duct  in  uniform 
flow  using  panels  set  at  varying  pitch  angles.  The  resulting 
chordwise  pressure  distribution  should,  of  course,  be  un¬ 
changed. 

Table  2  Effect  of  the  number  of  panda  on  the  computed  drag 
of  an  axisymmetric  duct  In  Imriecid  flow.  The  area  uaed  in 
defining  the  drag  coefficient  In  this  case  la  half  of  the  total 
wetted  surface  area  of  the  duct  The  exact  value  of  the  drag 
coefficient  la  zero 


N 

M 

Co 

Is 

9 

0  0121 

36 

IS 

0  0081 

36 

36 

00070 

60 

60 

0.0036 

Our  initial  results,  on  the  contrary,  showed  erratic  be¬ 
havior  and  poor  convergence  as  the  pitch  angle  of  the 
duct  paneling  was  decreased.  This  turned  out  to  be  due 
to  two  effects — the  first  being  the  fact  that  helical  panels 
become  very  nonplanar  as  the  pitch  is  decreased.  This  is 
very  evident  from  a  comparison  of  Fig.  13  and  Fig.  15. 
The  solution  to  this  part  of  the  problem  was  to  divide  the 
quadrilateral  panels  into  pairs  of  triangular  panels,  each 
with  the  same  uniform  source  and  dipole  strength.  The 
number  of  unknowns  was  therefore  the  same  as  with  quad¬ 
rilateral  panels. 

The  second  effect  was  the  sensitivity  of  the  potential  to 
the  radial  position  of  the  control  point  on  the  nonplanar 
panels.  We  found  that  the  best  position  could  be  found 
from  the  following  simple  argument.  Suppose  that  we 
wanted  to  compute  the  surface  potential  induced  by  a 


Fig.  15  Panel  arrangement  for  the  duct  shown  In  Fig.  13,  but  with 
helical  paneling  with  a  pitch  angle  of  20  deg 
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Fig.  16  Comparison  of  the  pressure  distributions  obtained  with  panels  oriented 
with  90-deg  pitch  as  shown  in  Rg.  1 3  and  with  panels  at  20-deg  pitch  as  shown 

in  Fig.  15 


uniform  distribution  of  normal  dipoles  located  on  the  seg¬ 
ment  of  a  cone  between  two  planes  normal  to  its  axis.  A 
discretized  version  of  this  problem  is  represented  by  any 
one  of  the  circumferential  strips  on  the  duct  surface  shown 
in  Fig.  13.  In  this  case  the  best  positions  for  the  control 
points  are  at  the  centroids  of  the  planar  quadrilaterals 
approximating  the  conic  surface.  In  particular,  the  radius 
of  a  control  point  is  less  than  the  radius  of  the  true  surface. 

But  we  also  know  that  the  dipole  panels  in  the  discre¬ 
tized  problem  are  equivalent  to  vortex  quadrilaterals,  and 
that  in  the  axisymmetric  case  the  meridional  elements 
cancel.  We  are  therefore  left  with  ring  vortices  at  each 
end,  which  in  the  discretized  problem  are  approximated 
by  polygons. 

Finally,  we  recognize  that  the  same  conclusion  holds  for 
helical  panels.  The  potential  induced  at  a  fixed  point  will 
be  unchanged  if  the  rings  are  connected  by  helical  lines 
of  arbitrary  pitch  and  the  problem  is  reformulated  in  terms 
of  dipole  panels.  Thus  we  conclude  that  the  correct  radial 
location  of  the  control  point  is  at  the  centroid  of  the  equiv¬ 
alent  planar  meridional  panel. 

The  numerical  justification  for  this  argument  may  be 
found  in  Fig.  16,  where  it  can  be  seen  that  the  pressure 
distributions  obtained  with  90-deg  and  20-deg  pitch  angles 
are  nearly  identical.  We  have  found  this  to  be  true  for 
several  different  duct  geometries  which  we  have  explored, 
whereas  the  results  obtained  with  other  choices  of  position 
of  the  control  points  were  not  nearly  as  satisfactory. 

Additional  geometrical  manipulations  are  required  to 
provide  compatibility  between  the  propeller  and  duct 
grids.  The  propeller  blade  tip,  and  its  extension  into  the 
trailing  vortex  wake  downstream,  is  adjusted  to  follow  the 
interior  contour  of  the  duct  with  a  prescribed  radial  gap. 
Finally,  the  axial  spacing  on  the  hub  and  duct  is  deter¬ 
mined  in  such  a  way  as  to  match  the  grid  on  the  blade. 
This  requires  cosine  spacing  on  the  interior  of  the  duct 
from  the  propeller  leading  edge  to  trailing  edge.  Two 


additional  regions  of  cosine  spacing  are  introduced,  ex¬ 
tending  from  the  duct  leading  edge  to  the  leading  edge 
of  the  propeller,  and  from  the  propeller  trailing  edge  to 
the  trailing  edge  of  the  duct.  Finally,  the  propeller  trailing 
vortex  wake  grid  is  adjusted  to  match  the  axial  spacing  on 
the  duct  and  on  the  duct  wake.  This  arrangement  is  similar 
to  that  developed  by  Van  Houten  [4]  for  use  with  his  vortex 
lattice  ducted  propeller  code. 

The  final  result  is  the  panel  extravaganza  depicted  in 
Fig.  17.  The  propeller  shown  has  five  blades,  and  a  nominal 
pitch  /  diameter  ratio  of  1 .4.  The  blade  outline  is  similar 
to  that  employed  in  the  Wageningen  KA  series  [32].  How¬ 
ever,  to  simplify  its  geometrical  description,  both  the  duct 
and  propeller  are  derived  from  a  single  thickness  form 
and  mean  line  with  maximum  thickness  and  camber 
matching  those  of  the  actual  KA  series.  Following  common 
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practice,  we  have  used  a  DTNSRDC  modified  NACA-66 
thickness  form  [33]  and  an  NACA  a  =  0.8  mean  line  [34]. 

Since  the  loading  on  the  duct  is  not  Circumferentially 
uniform,  the  wake  dipole  sheet  must  also  vary  in  strength. 
A  helical  pattern  for  the  duct  wake  is  assumed,  and  the 
pitch  of  this  wake  matches  the  tip  vortex  pitch  of  the 
propeller.  For  the  present,  we  have  suppressed  the  con¬ 
traction  of  the  propeller  tip  vortex.  The  duct  wake  pro¬ 
ceeds  downstream  with  a  constant  radius  equal  to  that  of 
the  duct  trailing  edge,  and  the  propeller  tip  vortex  pro¬ 
ceeds  downstream  with  the  same  radius  reduced  by  the 
tip  gap.  This  is  obviously  an  idealization  since  the  vortex 
systems  from  the  propeller  and  duct  must  actually  interact 
in  a  very  complicated  way.  However,  it  is  assumed  that 
this  wake  model  will  yield  reasonably  realistic  values  of 
induction  on  the  propeller  and  duct. 

For  the  steady  flow  problem,  the  loading  on  each  blade 
and  on  each  inter  blade  sector  of  the  hub  and  duct  is  the 
same.  Thus,  as  with  the  propeller,  the  number  of  un¬ 
knowns  is  equal  to  the  total  number  of  blade,  duct  and 
hub  panels  divided  by  the  number  of  blades.  In  addition, 
an  important  symmetry  exists  due  to  the  axisymmetric 
geometry  of  the  duct  and  hub.  If  the  angular  spacing  of 
the  panels  on  both  the  hub  and  duct  is  uniform  and  equal, 
the  influence  coefficients  will  depend  only  on  the  differ¬ 
ence  in  the  indices  identifying  the  helical  strip  containing 
the  panel  and  containing  the  control  point.  This  reduces 
the  number  of  influence  coefficients  to  be  computed  by 
a  factor  of  the  number  of  circumferential  panels  between 
blades.  Since  this  number  needs  to  be  on  the  order  of  10 
or  20,  the  reduction  in  computing  time  and  storage  is 
substantial. 

In  order  to  take  advantage  of  this  symmetry,  the  number 
of  circumferential  panels  on  the  duct  and  hub  must  be 
equal.  This  may  result  in  finer  spacing  than  necessary  on 
the  hub,  but  the  overall  reduction  in  computing  time  in¬ 
troduced  by  this  symmetry  far  outweighs  the  time  which 


would  be  saved  by  removing  a  few  circumferential  panels 
on  the  hub. 

As  with  the  propeller  panel  code,  an  accelerated  iter¬ 
ative  matrix  solver  is  used  for  the  duct  and  hub  solution. 
The  velocity  induced  by  the  duct  and  hub  at  the  propeller 
control  points  is  then  computed  directly  from  the  panel 
influence  functions.  This  approach  is  accurate  since  most 
of  the  field  points  are  far  from  the  duct  and  hub  in  terms 
of  panel  dimensions,  and  the  nearest  po  nts  are  located  on 
a  properly  aligned  grid. 


7.  Ducted  propeller  results 

We  will  next  discuss  some  of  the  results  obtained  for  the 
ducted  propeller  illustrated  in  Fig.  17.  Following  the  con¬ 
clusions  of  Appendix  4,  we  have  made  calculations  for  two 
extreme  limits  of  the  tip  gap — a  gap  of  4  percent  of  the 
radius,  for  which  a  potential  flow  model  of  the  tip  and  the 
hub  is  reasonably  valid,  and  zero  gap,  where  no  gap  flow 
exists. 

Figure  18  shows  how  the  iterations  between  the  duct 
and  blade  solution  converge  in  the  case  of  the  4  percent 
gap.  The  first  estimate  of  the  radial  distribution  of  circu¬ 
lation  on  the  blades  is  high,  since  at  this  stage  none  of  the 
propeller  circulation  has  been  transferred  to  the  duct.  In 
the  second  iteration,  the  duct  circulation  induced  by  the 
propeller  has  resulted  in  increased  flow  speed  through  the 
duct,  thus  decreasing  the  loading  on  the  blades.  This  effect 
is  slightly  overestimated,  since  the  first  iteration  blade 
loading  was  too  high.  By  the  third  iteration,  the  results 
have  practically  converged,  and  almost  no  difference  can 
be  seen  in  the  figure  between  the  fourth  and  seventh 
iterations.  The  same  pattern  exists  for  the  solution  on  the 
duct. 

Some  indications  of  sensitivity  to  numbers  of  panels  is 
given  in  Figs.  19  and  20.  The  first  shows  the  circumfer- 
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Fig.  IS  Illustration  of  the  iterative  procedure  used  to  obtain  the  duct  and  blade 
solutions.  The  circulation  distribution  on  the  blade  obtained  for  each  of  seven 
iterations  is  shown 


12 


Ducted  Propellers 


POSITION  BETWEEN  BLADES-DEG 

Fig.  19  Example  of  the  influence  of  numbers  of  panels  on  the  calculated  cir¬ 
culation  on  the  duct 


ential  variation  in  total  circulation  on  the  duct  obtained 
with  different  numbers  of  chordwise  and  circumferential 
panels.  The  two  finest  spacings  produce  nearly  identical 
results,  while  the  coarsest  spacing  yields  results  which  are 
about  3  percent  low.  Figure  20  shows  the  chordwise  pres¬ 
sure  distribution  on  the  duct  for  one  of  the  panels  closest 
to  the  blade  tip.  The  differences  in  pressure  distribution 
obtained  with  60,  80  and  100  panels  is  due  principally  to 
the  fact  that  the  actual  position  on  the  duct  moves  closer 
to  the  blade  tip  as  the  number  of  panels  is  increased. 


Figure  21  shows  the  radial  distribution  of  circulation  on 
the  blades  for  zero  and  4  percent  gap,  and  the  trends  near 
the  blade  tip  are  as  expected.  The  fact  that  the  maximum 
circulation  increases  as  the  gap  is  increased  is  evidently  a 
result  of  the  coupling  of  the  blade  and  duct  circulation 
altering  the  mean  flow  speed  through  the  duct.  Figure  22 
shows  the  effect  of  advance  coefficient  on  the  radial  dis¬ 
tribution  of  circulation  on  the  blades  for  the  case  of  zero 

gap- 

Figure  23  shows  the  effect  of  the  gap  on  the  circum- 


Flg.  20  Example  of  the  influence  of  numbers  of  panels  on  the  calculated  pressure 
distribution  on  the  duct 
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Fig.  21  Influence  of  tip  gap  on  the  radial  distribution  of  circulation  on  the  blades 


ferential  distribution  of  circulation  on  the  duct.  Decreasing  circulation  on  the  blades  at  a  sequence  of  positions  over 
the  gap  from  four  percent  to  zero  greatly  increases  the  the  radius.  The  blades  are  evidently  operating  with  con- 
mean  circulation,  as  well  as  its  circumferential  variation,  siderable  angle  of  attack  loading  in  this  condition. 

In  addition,  the  zero  gap  case  shows  a  discontinuity  in  Finally,  Figs.  25  and  26  show  the  chordwise  pressure 
circulation  at  the  blade  tip,  which  is  reasonably  close  to  distributions  on  the  duct  for  each  of  the  helical  strips  of 
the  value  of  the  finite  circulation  at  the  blade  tip.  panels  between  a  pair  of  blades.  In  the  case  of  the  4  percent 

We  turn  next  to  detailed  flow  distributions  on  the  blade  gap,  the  pressure  minimum  is  near  the  blade  tip,  but  is 
and  duct.  Figure  24  shows  the  choidwise  distribution  of  continuous  in  passing  across  the  blade  tip.  On  the  other 
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circulation  distributions  on  the  biade  at  an  advance  coefficient  Fig.  26  Pressure  distribution  on  the  interior  surface  of  the  duct  for  the  case  of 

j  =  o.e  with  zero  tip  gap  J  =  0.7  and  zero  tip  gap.  The  coordinates  are  the  same  as  in  Fig.  25 


Fig.  27  Pressure  distribution  on  the  interior  surface  of  the  duct  for 
the  case  of  J  —  0.6  and  4  percent  tip  gap.  The  coordinates  are  the 
same  as  in  Fig.  26 


hand,  with  zero  gap,  the  pressure  minimum  on  the  duct 
is  on  the  suction  side  of  the  blade  tip,  and  the  pressure 
coefficient  is  discontinuous  across  the  tip.  This  disconti¬ 
nuity  is,  of  course,  related  to  blade  loading  as  evidenced 
by  the  difference  between  Figs.  27  and  26.  The  calculated 
total  propeller  and  duct  forces  for  these  cases  are  given 
in  Table  3. 
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Fig.  28  Notation  for  a  general  body  for  the  application  of  Green’s 
theorem 


Review  of  the  basic  theory 

The  basic  mathematical  theory  behind  the  various  panel  meth¬ 
ods  will  be  summarized  and  the  characteristics  of  each  method 
will  be  compared  in  order  to  determine  the  most  suitable  method 
for  the  analysis  of  marine  propellers  and  ducts.  The  common  basis 
of  the  apparently  different  panel  methods  will  then  become  ev¬ 
ident. 

Panel  methods  can  be  grouped  as  follows: 

•  Potential  field  formulation: 

— perturbation  potential  method 
— total  potential  method. 

•  Velocity  field  formulation: 

— mixed  source  and  dipole  method 
— dipole  method  (equivalently,  vorticity  method) 

— source  based  method 


•  A  Kutta  condition  is  required  at  the  trailing  edge  to  uniquely 
specify  the  circulation.  In  its  most  general  form,  it  states  that  the 
flow  velocity  at  the  trailing  edge  remains  bounded,  that  is 

<  <*>  (9) 

•  On  the  outer  control  surface  S„,  the  perturbation  velocity  due 
to  the  body  should  vanish  in  the  limit  where  this  surface  is  an 
infinite  distance  from  the  body. 

V<f>  0,  as  S„  -*  oo  (10) 

According  to  Lamb  [35],  this  boundary-value  problem  from  the 
velocity  potential  outside  the  body  surface  can  be  transformed 
into  an  integral  equation,  upon  consideration  of  a  fictitious  fluid 
in  V,  which  is  the  domain  internal  to  the  body  surface  S,: 


Statement  of  the  problem 

Consider  a  closed  three-dimensional  domain  V  with  boundary 
S,  the  unit  normal  vector  n  to  S  being  oriented  into  V,  as  shown 
in  Fig.  28.  The  boundary  5  is  composed  of  the  body  surface  S„ 
the  wake  surface  Sw,  and  the  outer  control  surface  S.  surrounding 
the  body  and  wake  surface.  The  body  is  subject  to  the  inflow 
velocity  {/_.  With  the  assumption  that  the  fluid  in  V  is  incom¬ 
pressible,  inviscid  and  irrotational,  there  exists  a  perturbation 
velocity  potential  <fr,  which  satisfies  the  Laplace  equation 

V*<|(  =  0  (4) 


fa 

_  /SMq)  _  d<l»'(<?)\  1  1  ds 

V  dn,  dn,  )  ff(p#)J 

+  f  f  &4>(q)  JL —1— dS 
J  J  Hn,R(pw) 

Sw 


(11) 


A  boundary-value  problem  can  be  constructed  by  specifying 
boundary  conditions  on  the  boundary  S  as  follows: 

•  The  kinematic  boundary  condition  should  be  satisfied  on  the 
solid  body  surface  S,; 

g  n  (5) 

on 

•  The  wake  surface  Sw  is  assumed  to  have  zero  thickness.  The 
normal  velocity  jump  and  the  pressure  jump  across  Sw  is  zero, 
while  a  jump  in  the  potential  is  allowed: 


where 


4  =  perturbation  velocity  potential  in  V 
4'  =  perturbation  velocity  potential  in  V 
p(x,y,z)  =  field  point  where  induced  potential  is  calculated 
<?((,T),{)  =  source  point  where  singularity  is  located 
R(pw)  =  distance  between  points  p  and  q 


d_ 

9m 


V<*  -  €>*  +  (If  -  11)*  +  (*  -  0* 


normal  derivative  with  respective  to  point  q 


<ApLjw  -  P*  -  P~  “  0  (6) 


For  the  steady  lifting  problem,  the  potential  jump  across  the  wake 
surface  is  the  same  as  the  circulation  around  the  body,  and  is 
constant  in  the  stream  wise  direction  on  Sw: 

<*♦)-.„  T  (8) 


This  equation  may  be  regarded  as  a  representation  of  the  ve¬ 
locity  potential  in  terms  of  a  normal  dipole  distribution  of  strength 
(4  -  ♦')  on  the  body  surface  S*  a  source  distribution  of  strength 
(94l9n)  -  (94'  I9n)  on  S*,  and  a  normal  dipole  distribution  of 
strength  44  on  the  wake  surface  Sw. 

Because  the  fictitious  fluid  inside  S,  does  not  have  physical 
meaning,  we  can  choose  the  internal  velocity  potential  4'  to  suit 
our  convenience.  By  choosing  a  different  4'  in  equation  (11),  we 
can  formulate  a  different  kind  of  the  panel  methods  which  uses 
a  different  set  of  the  singularities. 
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Potential  field  formulation 

If  we  choose  the  fictitious  potential  as  <£'  =  0  on  S„  equation 
(11)  for  the  field  point  p  on  the  body  surface  SB  becomes 


2*4*p)  =  JJ[*<„)A_±_ 

-1—1 
(p;<?)l 


dn,  R(p;q) J 


(12) 


SS 


A<f>(<7> 


1 


dn,  R(p;q) 


dS 


Here  the  surface  integral  on  Sa  must  be  defined  to  exclude  the 
immediate  vicinity  of  the  singular  point. 

Because  ctyldn  is  known  on  SB  from  the  boundary  condition, 
equation  (5),  equation  (12)  is  a  Fredholm  integral  equation  of  the 
second  kind  for  the  dipole  strength  4>,  which  is  also  the  potential 
value  on  the  body  surface  SB.  The  potential  jump  across  the  wake 
surface  can  be  set  to  the  difference  between  the  potential  values 
of  the  upper  and  lower  surface  at  the  trailing  edge,  which  replaces 
the  Kutta  condition.  Discretization  of  equation  (12)  will  lead  to 
a  linear  system  of  equations  for  the  unknown  4>.  The  surface 
velocity,  hence  the  pressure,  on  SB  can  be  calculated  by  a  nu¬ 
merical  differentiation  of  the  potential  distribution.  This  form  of 
a  panel  method  was  first  used  by  Morino  (16],  and  adopted  in 
the  present  paper.  We  will  refer  to  this  as  Morino’s  method  or 
the  perturbation  potential  method. 

_  If  there  exists  the  inflow  velocity  potential,  such  that  7<#>_  = 
[/„,  we  can  formulate  another  form  of  the  panel  method  by 
choosing  the  internal  potential  in  equation  (11)  as  a  negative  of 
the  inflow  velocity  potential,  that  is,  =  —  <b„  The  source 
strength  in  equation  (11)  becomes  zero  because  of  the  boundary 
condition,  equation  (5),  while  the  dipole  strength,  which  is  the 
difference  between  <j>  and  <f>',  becomes  the  total  potential 


<t>  —  <j>—  <f>  +  4,«='!> 


(13) 


As  the  point  p  approaches  the  body  surface  S„  the  contribution 
from  the  immediate  surface  S,  on  SB  in  the  first  term  of  equation 
(11)  is 


lim  lim  f  f  (4,  -<(,')  J-  — L_  dS 
s  -or-sB  J  J  dn R(p;q) 


~  2ir($  —  <£')  =  2 »r4>  (14) 


The  resulting  equation  is 

2ir<I>(p)  =  4ir4>.(p)  + 


SS 


*(</) 


1 


dn,  R{p;q) 


dS 


SI 


(15) 


A4>(<?) 


dn,  R(PW) 


dS 


This  equation  can  be  regarded  as  a  representation  of  the  total 
velocity  potential  in  terms  of  a  normal  dipole  distribution  only 
on  the  body  surface  S,  and  the  wake  surface  Sw.  Given  the  inflow 
velocity  potential  values,  this  is  also  a  Fredholm  integral  equation 
of  the  second  kind  for  the  total  potential  <t>  Discretization  of  this 
equation  gives  another  form  of  panel  method,  which  we  will  refer 
to  as  the  total  potential  method. 


Velocity  field  formulation 

Instead  of  forming  an  integral  equation  in  the  potential  field, 
we  can  alternatively  construct  one  in  the  velocity  field.  Take  a 
normal  derivative  of  equation  (11)  with  respect  to  the  field  point 
p;  then  the  resulting  equation  when  the  field  point  p  is  on  S,  is 


4*  +  /  /  s;  fas) iS 


ss 


n(q) 


l 


dnjdn,  R(p;q) 


dS 


(16) 


-ss 


a<m<?) 


3s 


dn,dn,  R(p;q) 


dS 


where 


dn  dn 


and  p  =  <t>  —  <f>' 


Because  the  velocity  induced  by  the  dipole  distribution  is  equiv¬ 
alent  to  that  induced_by  a  vorticity  distribution  on  the  same 
surface  with  strength  y ,  which  is  calculated  as  a  vector  product 
of  the  normal  vector  and  the  local  surface  gradient  of  the  dipole 
strength  [11],  we  can  write  alternatively: 


4*  ^£7  =  +  !  I  £  fir) ds 

SB 

+  |  J  n,-Y(<?)  X  7p(~)rf5  (17) 

+  J  J  np-y(q)  X 

sw 

where 

y  =  n,  X  Vs^(4>  ~  4>') 

Here  again,  as  we  choose  a  different  value  for  the  internal 
potential  4>’  inside  S„,  we  can  express  the  normal  velocity  on  S„ 
in  terms  of  a  different  set  of  singularities. 

If  we  choose  the  internal  potential  in  equation  (16)  as  4>'  =  0, 
it  can  be  shown  that  34»'  /  dn,  =  0  on  SB;  then 

+  J7'l(<,)*ksdS  (18) 

SB 

+  J7 

sw 

where 


<r  =  --  and  p  *=  $ 


This  can  be  regarded  as  the  normal  induced  velocity  at  p  due 
to  the  mixed  distribution  of  normal  dipoles  of  strength  <j>,  the 
source  distribution  of  strength  d$  I  dn  on  SB  and  the  normal  dipole 
distribution  of  strength  A4>  on  Sw.  Given  the  source  strength  and 
the  left  hand  side  of  the  equation  from  the  body  boundary  con¬ 
dition,  this  is  an  integral  equation  of  the  first  kind  for  the  unknown 
dipole  strength  p.  Discretization  of  this  equation  will  give  another 
form  of  panel  method,  which  we  will  refer  to  as  the  mixed  source 
and  dipole  method.  This  can  be  regarded  as  the  velocity  field 
formulation  of  the  perturbation  potential  method. 

An  equivalent  formulation  derived  from  equation  (17)  leads  to 
a  mixed  distribution  of  sources  and  vortices  instead  of  dipoles 
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2ir 


SB 

+  j  j  np-y(Q)  X  7„  (jpj  dS 

SB 

+  J  /  *<<7)  X  7”  (V)  dS 


(19) 


tailed  investigation  that  led  to  this  decision  is  described  in  [17]. 
In  this  Appendix  we  will  describe  the  implementation  of  the 
pressure  Kutta  condition  on  a  duct  operating  in  conjunction  with 
a  propeller.  The  implementation  of  the  Kutta  condition  in  the 
present  propeller  panel  code  is  very  similar. 

Green’s  formula  (11)  for  the  perturbation  potential  on  the  duct 
is  discretized  as  follows: 


I  EMb  =  I  Su[U£.nt]  +  <*>'+  X  W_,  • 


(23) 


where 

dd>  -  _ 

cr  =  —  and  y  =  n,  X  Vi-A 
an 

If  we  choose  the  internal  potential  as  3<j>'  I  bn  =  ckj>  /  3n,  the 
source  strength  becomes  zero  and  equation  (16)  becomes 


4ir 


*)>(P) 

dn„ 


”  J  J  ^ 


dn.dn,  ft 


idS 


J  l**®  si  7s  » 


where  p.  =  <(>  +  <j>.  =  4>. 

This  is  an  integral  equation  of  the  first  land  for  the  unknown 
dipole  strength  p.  The  panel  method  derived  from  equation  (20) 
will  be  referred  to  as  the  dipole  method.  This  can  be  regarded 
as  the  velocity  field  formulation  of  the  total  potential  method. 

Due  to  the  equivalence  between  dipoles  and  the  vortices,  equa¬ 
tion  (20)  can  be  written  in  a  different  form. 


Air 


&HP) 


y(q)  X  7, 


(t) 


dS 


/;*■ 


y(q)  x  7, 


% 


where  D„  and  are  the  potentials  at  the  control  point  i  due  to 
unit  constant  strength  dipoles  and  sources,  respectively,  placed 
at  panel  j,  and  <f>'  is  the  perturbation  potential  due  to  the  pro¬ 
peller.  The  quantity  A$„  represents  the  jump  in  the  potential  at 
the  m*  strip  on  the  duct  wake,  and  Wm  is  the  influence  of  the 
mth  dipole  wake  strip  at  the  i*  control  point. 

The  unknown  &4>„  will  be  determined  by  applying  the  Kutta 
condition  at  the  duct  trailing  edge. 

The  Kutta  condition  requires  that  the  velocity  at  the  trailing 
edge  of  the  duct  be  finite.  In  the  numerical  formulation  of  the 
problem,  we  will  implement  the  Kutta  condition  by  requiring 
that  the  pressures  at  the  upper  and  lower  control  points  at  the 
trailing  edge  be  equal.  This  can  be  expressed  as  follows: 


Ap„  =  p„ 


0 


for  m  =  1,  Af* 


(24) 


dS  (21) 


where  y  =  J,  X  72  ^<t>.  We  will  refer  to  the  panel  method  from 
this  equation  as  the  vorticity  method. 

If  we  choose  the  vorticity  strength  in  equation  (17),  such  that 
it  has  a  given  shape  function  g(t)  along  the  chordwise  panels,  and 
the  spanwise  circulation  T(j)  is  yet  to  be  determined,  then 

4w  ^7  =  2™{p)  +  J  J  *(q)  £r,  Gr) ds 

SB 

+JJvyx7  (22) 

hi 

+J7  x7'®</s 


where  y  =  T(s)g(f)  t  y,  t  y  is  the  direction  of  the  vorticity,  and 
s  and  t  are  the  spanwise  and  chordwise  coordinates.  Given  the 
normal  velocity  on  S,  from  the  boundary  condition  and  the  vor¬ 
ticity  shape  function  g(t),  this  is  a  Fredholm  integral  equation  of 
the  second  kind  for  the  unknown  source  strength  t t  and  the 
spanwise  circulation  distribution  T(s).  This  is  the  form  of  the 
original  surface  source  method  by  Hess  [31].  We  will  refer  to  this 
form  as  the  source  based  method. 


Appendix  2 

Implementation  of  the  pressure  Kutta  condition 

As  stated  in  the  beginning  of  the  paper,  the  present  panel 
method  employs  an  explicit  pressure  Kutta  condition.  The  de- 


where  M  is  the  total  number  of  circumferential  panels  on  the 
duct,  and  A/*  is  the  total  number  of  circumferential  panels  be¬ 
tween  two  blades. 

A  direct  solution  of  the  resulting  system  of  equations  (23)  and 
(24)  is  difficult  due  to  the  nonlinear  character  of  equation  (24). 
Therefore  an  iterative  solution  algorithm  is  employed.  At  the  k"' 
iteration,  we  solve  the  linear  system  of  equations  (23)  with  the 
values  of  A$  J*1  determined  from  the  (k  —  1)”  iteration.  The  values 
of  ApJ*1  are  given  by  (24),  with  the  values  of  the  pressures  pmv 
and  pmL  determined  as  described  in  Appendix  3.  If  ApJ*1  is  not 
equal  to  zero  within  the  desired  tolerance,  we  proceed  to  another 
iteration  with  A«t>J4^"  determined  as  follows: 

[A*]'**11  =  [A<f.r  -  m  1  •  [Ap]'*'  (25) 

[Ap]  =  [Ap,,  Ap, . A pA/»]r  (26) 

[A<H  =  [A*„  A<f>, . Ad >M„]T  (27) 


where 


and  [ J ]  1  is  the  inverse  of  the  Jacobian  matrix,  the  elements  of 
which  are  defined  as 


/  -  d<AP‘> 

"  d(A  <M 


(28) 


with  the  values  of  the  partial  derivatives  approximated  numeri¬ 
cally  as 


<KA p,)  Ap,1*'  -  Ap,'0 


*A<M  ~  A<f>;*‘  -  A<f>;° 


(29) 


where  A p,10’  corresponds  to  the  initial  guess  A<f>,10’,  and  Ap,1*1  cor¬ 
responds  to  A  sty*’,  a  perturbation  to  the  initial  guess  defined  as 


and 


A*;*’  =  (1  -  d)A<f.'° 


A ♦/*’  =  Ad>,'0'  fox  l*  j 


(30) 


(31) 


where  0  is  a  small  number.  We  have  selected  0  =  0.01,  although 
the  algorithm  has  been  found  to  converge  for  a  large  range  of  fi. 

The  initial  guess,  A^J01,  is  taken  as  the  difference  of  the  po¬ 
tentials  at  the  upper  and  lower  control  points  at  the  trailing  edge 
of  the  duct 


A*J01  = 


(32) 


The  initial  guess  is  therefore  the  original  Morino  [16]  Kutta  con¬ 
dition. 
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Fig.  29  Comparison  of  duct  chordwise  pressure  distributions  obtained  before 
and  after  the  application  of  the  iterative  pressure  Kutta  condition 


Figure  29  shows  a  typical  chordwise  pressure  distribution  on 
a  duct  before  and  after  the  application  of  the  iterative  Kutta 
condition.  The  chordwise  panel  shown  in  this  example  is  adjacent 
to  the  blade  tip,  where  crossflow  effects  are  greatest.  It  is  evident 
that  the  pressures  at  the  trailing  edge  are  initially  unequal,  but 
become  equal  after  application  of  the  iterative  Kutta  condition. 
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where  Ut  is  the  relative  effective  incoming  velocity  and  7<J>  is 
the  perturbation  velocity  at  the  point  under  consideration. 

To  determine  the  pressure  distribution  on  the  body,  it  is  nec¬ 
essary  to  determine  q,  at  the  control  points  on  its  surface  If  n 
is  the  normal  unit  vector  on  the  surface  of  the  body,  then  ac¬ 
cording  to  the  kinematic  boundary  condition  on  the  bod) 

q,  ■  n  —  0  (37) 

By  combining  equations  (36)  and  (37),  we  obtain  the  result 


Computation  of  pressure  distributions 

After  the  discrete  values  of  the  potential  on  the  body  have 
been  computed  the  pressure  distribution  on  the  surface  can  be 
determined  by  applying  Bernoulli's  equation.  Wg  assume  that  the 
body  rotates  with  a  constant  angular  velocity  n  since  this  has  a 
direct  application  to  the  propeller  or  to  the  duct  problem.  Calling 
p  the  pressure,  q,  the  total  fluid  velocity^  with  respect  to  the 
rotating  system,  p  the  fluid  density,  and  7  the  vector  distance 
from  the  axis  of  rotation,  we  define  the  total  pressure  head,  H, 
according  to  Batchelor  [36]  as  follows 

H  X  7,1  <33) 

The  flow  field  around  a  ducted  propeller  in  axisymmetric  flow 
will  be  steady  when  referred  to  a  coordinate  system  rotating  with 
the  propeller.  In  that  case  Bernoulli’s  equation,  as  shown  in  Batch¬ 
elor  [36],  can  be  expressed  as 


7  H  =  q,x(w  +  2ft)  (34) 

where  <o  is  the  vorticity  of  the  fluid  with  respect  to  the  rotating 
system.  The  term  3+20  represents  the  vorticity  vector  with 
respect  to  the  absolute  system.  Equation  (34)  requires  that  the 
total  pressure  head  along  stream  surfaces,  Hs,  remain  constant: 


Hs  =  constant  =  —  +  i  l/.* 
P  2 


(35) 


where  U.  is  the  axial  inflow  velocity  and  p.  the  pressure  far 
upstream  on  the  same  stream  surface. 

In  the  case  of  uniform  inflow,  H  will  be  constant  everywhere 
in  the  flow  field. 

The  total  velocity  vector  q ,  can  be  decomposed  as  follows: 


Q,  =  U*  + 


(36) 


where 


q,  =  Vt  -  ( ti  ■  V, )  n  +  7s<k 
7s<f  =  7<J>  -  (  n  -  7*)K 


(38 1 
(39) 


The  term  7<4>  corresponds  to  the  component  of  the  perturba¬ 
tion  velocity  tangent  to  body  surface.  If  we  call  7  and  I  the 
unit  vectors  tangent  to  the  surface  along  the  two  grid  directions, 
the  following  equations  hold: 


cW> 

ds 

a* 

at 


*  •  7,<f 

(40) 

7  -  7 

(41) 

By  using  equations  (40)  and  (41)  we  can  express  7,<J>  as: 


7S4>  = 


a<t>  -  _  — .  as  - 

II*  x  0|s 


-(*•/)*] 


(42) 


The  quantities  d^/ds  and  d$!dt  correspond  to  the  projections 
of  the  perturbation  velocity  along  the  two  grid  directions,  and 
these  will  be  obtained  numerically. 

In  the  special  case  of  a  duct,  d$ids  is  determined  from  the 
second-order  finite  difference  of  the  computed  values  of  the  po¬ 
tentials  at  the  control  points  along  each  chordwise  duct  station. 
d$ldt  is  determined  from  the  derivative  of  the  cubic  spline  pass¬ 
ing  through  the  potentials  along  each  circumferential  duct  station. 

The  values  of  thepressures  are  computed  by  using  equations 
(35),  (38)  and  (39).  The  pressures  are  finally  expressed  in  terms 
of  the  pressure  coefficient  C,,  which  is  defined  as: 


Cr~ 


(43) 
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Appendix  4 


Modeling  of  the  tip  gap 

The  radial  clearance  between  the  propeller  and  the  interior 
surface  of  the  duct  is  generally  small  compared  with  the  overall 
dimensions  of  the  device.  As  a  result,  a  large  mismatch  in  length 
scale  exists  between  the  overall  flow  around  the  duct,  hub,  and 
propeller  blades  and  the  local  flow  through  the  tip  gap.  This  would 
not  be  a  serious  problem  except  for  the  fact  that  the  local  gap 
flow  has  a  major  influence  on  the  global  flow.  A  similar  problem 
is  encountered  in  the  analysis  of  axial  flow  turbomachines. 

The  flow  around  an  isolated  duct,  hub  or  set  of  propeller  blades 
can  be  represented  with  reasonable  accuracy  as  a  potential  flow. 
What  really  happens  when  the  tip  of  a  propeller  blade  is  brought 
in  close  proximity  to  the  duct  is  not  obvious.  The  characteristics 
of  a  potential  flow  analysis  can  be  seen  most  easily  by  examining 
the  solution  for  a  simple  lifting  line  of  span  s  and  constant  unit 
downwash  u>*  =  1  operating  at  a  span  wise  distance  h  from  an 
infinite  plane  wall.  The  wall  can,  of  course,  be  represented  by  an 
image  lifting  line.  The  span  wise  distribution  of  circulation 


Git)  = 

tv’s 


and  the  distribution  of  velocity  on  the  wall  and  in  the  gap  region 
can  be  computed  using  a  vortex  lattice.  The  number  of  vortex 
lattice  panels  required  for  accurate  results  increases  as  the  gap 
ratio  h/s  decreases,  but  it  is  easy  to  obtain  graphically  exact  results 
with  very  little  computing  effort. 

Figure  30  shows  the  spanwise  distribution  of  circulation  ob¬ 
tained  for  various  gap  ratios.  The  lowest  curve  is  for  a  relatively 
large  gap  ratio  of  0.04,  while  the  next  four  curves  show  how  the 
circulation  changes  as  the  gap  is  successively  halved.  Finally,  the 
top  curve  represents  zero  gap,  where  the  effect  of  the  wall  image 
is  to  create  the  same  flow  as  for  an  elliptically  loaded  lifting  line 
of  span  2s.  Presenting  the  results  in  this  way  emphasizes  the  fact 
that  a  tiffing  line  with  an  e  tremely  small  gap  does  not  behave 
at  all  like  one  with  zero  gap,  and  that  an  infinitesimal  change  in 
the  gap  produces  a  large  global  change  in  the  circulation  distri¬ 
bution. 

The  reason  for  this  can  be  found  in  Fig.  31,  which  shows  the 
distribution  of  velocity  along  the  wall  induced  by  the  lifting  line 
and  it's  image.  Here  we  see  that  the  lifting  line  with  a  gap  ratio 
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Rfl.  31  Nondimensional  velocity  w/w'  along  the  wall  induced  by 
a  lifting  line  and  its  image  for  the  same  gap  ratios 


of  0.0025  induces  a  peak  velocity  of  over  50  times  the  downwash 
velocity!  While  the  gap  is  negligibly  small,  the  existence  of  such 
hurricane  force  flow  speeds  explains  why  the  circulation  is  af¬ 
fected  everywhere. 

Figure  32  shows  the  volumetric  flow  through  the  gap  as  a 
function  of  gap  ratio.  While  the  gap  flow  eventually  goes  to  zero 
as  the  gap  is  closed,  it  is  approximately  constant  in  the  range  of 
practical  gap  ratios,  thus  indicating  that  the  gap  velocities  will  be 
roughly  inversely  proportional  to  the  gap  ratio. 

These  results,  however,  are  contrary  to  experimental  evidence, 
where  it  has  been  found  that  the  circulation  distribution  tends 
towards  the  zero  gap  limit  much  more  rapidly.  A  physical  expla¬ 
nation  is  that  the  cross  flow  around  the  tip  separates,  thus  forming 


a. 


1*10. 30  Nondktwngkxwl  spanwtse  circulation  distributions.  r/Wt,  Fig.  32  Nondiroenatonal  volumetric  flow/unit  chord  O/w’s  through 

for  a  lifting  line  next  to  a  wax  for  gap  ratios  h/s  =  0, 0.0025, 0.005,  the  gap  between  a  lifting  line  and  Its  image  as  a  function  of  the  g ap 
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a  leakage  vortex  sheet  shown  schematically  in  Fig.  33,  taken  from 
Lakshminarayana  [6],  This  is  fundamentally  similar  to  the  mech¬ 
anism  of  tip  vortex  generation  for  a  wing  or  hydrofoil  in  an  un¬ 
bounded  fluid.  However,  in  this  case  the  cross  flow  is  greatly 
increased  due  to  the  presence  of  the  wall,  and  the  flow  is  further 
complicated  by  the  presence  of  the  wall  boundary  layer. 

A  simple  model  of  the  crossflow  through  the  gap  which  appears 
in  several  sources  in  the  turbomachinery  literature  has  been  re¬ 
cently  reviewed  by  Van  Houten  [4].  In  this  model,  the  ieakage 
vortex  sheets  leaving  the  tip  of  the  physical  blade  and  its  image 
appear  as  the  boundaries  of  a  two  dimensional  jet  when  viewed 
in  a  plane  normal  to  the  blade  surface.  Thus  the  flow  is  approx¬ 
imated  by  one  through  an  orifice.  The  relationship  between  the 
volumetric  crossflow,  Q,  and  the  pressure  difference  driving  the 
flow  is 

Q=hCQ^2  <44) 

where  C0  is  a  discharge  coefficient  which  depends  principally  on 
the  tip  cross-sectional  shape,  the  wall  motion,  and  wall  boundary 
layer  characteristics.  Typical  values  compiled  by  Van  Houten  [4] 
range  from  0.76  to  0.92. 

The  recent  development  of  efficient  finite-difference  methods 
for  solving  the  two-dimensional  Navier-Stokes  equations  have 
made  it  possible  to  obtain  theoretical  crossflow  solutions  for  dif¬ 
ferent  tip  geometries  [7].  Such  calculations  provide  a  rational  basis 
for  determining  Cq,  but  do  not  alter  the  fundamental  nature  of 
the  basic  orifice  equation. 

One  can  make  an  order-of-magnitude  comparison  of  the  vol¬ 
umetric  flow  given  by  equation  (44)  and  the  lifting  line  results 
plotted  in  Fig.  32  by  equating  the  pressure  difference  Ap  in 
equation  (44)  to  the  lift  on  the  lift  per  unit  area  on  the  hydrofoil 
represented  by  a  lifting  line.  This  requires  the  introduction  of  an 
aspect  ratio  A  and  a  ratio  of  downwash  velocity  to  inflow  speed, 


Fig.  33  Schematic  of  gap  vortex  formation  around  a  blade  tip  and 
its  image,  from  Lakshminarayana  [6] 


w*/U.  The  orifice  equation  then  becomes 

(45> 

By  substituting  any  reasonable  values  for  the  circulation  and 
aspect  ratio  in  equation  (45),  one  concludes  that  for  small  gaps, 
the  mass  flow  driven  by  the  pressure  difference  across  the  gap  is 
orders  of  magnitude  smaller  than  the  potential  flow  value  derived 
from  lifting  line  theory.  On  the  other  hand,  for  large  gaps,  where 
the  idealization  of  two-dimensional  orifice  flow  becomes  less  valid, 
the  flow  predicted  from  (45)  is  of  the  same  order  of  magnitude 
as  the  lifting  line  value. 

While  the  ducted  propeller  flow  problem  is  more  complicated, 
the  trend  must  be  similar.  A  potential  flow  representation  of  the 
blades  and  duct  is  suitable  for  either  large  clearances  or  zero 
clearance,  where  there  is  no  flow  around  the  tip.  However,  for 
small  clearances  one  must  either  turn  to  a  complete  three-di¬ 
mensional  viscous  flow  solution  in  the  tip  gap  region,  or  adopt 
some  means  of  coupling  the  potential  flow  solution  to  a  local 
representation  of  the  viscous  crossflow  in  the  gap. 
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